Load required packages
# load pacman package to load or install other requried libraries
if (!require('pacman')) install.packages('pacman'); library(pacman)
# load (install if required) packages from CRAN
p_load("here","tidyverse","lubridate","janitor","plotly","doParallel")
# load/install packages from GitHub
p_load_gh("reichlab/zoltr", "reichlab/covidHubUtils")
Utilizing zoltr package to enable access to Zoltar API to the forecast archive and covidHubUtils package to that provide functions to read, plot and score forecast data.
Â
Importing observed data
We will use the Covid-19 cases as updated by CDC. We will download data from the data Table for daily Case Trends for The United States found here.
Also, we will pull daily COVID-19 cases from CDC at a state-level found here.
Data extracted September, 06, 2021. Here is the top rows of the incident COVID-19 cases data.
# csv export of daily cases from CDC COVID data tracker at national level
national_cases_daily <- read_csv(here("data", "data_table_for_daily_case_trends__the_united_states.csv"),
col_types = cols(Date = col_date(format = "%b %d %Y")),
skip = 2)
# csv export of daily cases from CDC COVID data tracker at state level
states_cases_daily <- read_csv(here("data", "United_States_COVID-19_Cases_and_Deaths_by_State_over_Time.csv"),
col_types = cols(submission_date = col_date(format = "%m/%d/%Y")))
head(national_cases_daily)
Â
Preparing data
We’ll aggregate the daily case counts into weeks, starting on Sunday as is the same as an epidemiological week (referred to as MMWR week, more info found here)
# clean up column names
national_cases_daily <- national_cases_daily %>%
clean_names()
# aggregate daily COVID cases into weekly case counts,
# start week on Sunday according to epidemiological week (MMWR week)
week <- as.Date(cut(national_cases_daily$date, "week", start.on.monday = FALSE))
national_cases_weekly <- aggregate(new_cases ~ week, national_cases_daily, sum)
head(national_cases_weekly)
# filter dates after August
national_cases_weekly <- national_cases_weekly %>%
filter(week < "2021-09-01")
The forecast models are pulled from the Zoltar forecast model archive using the covidHubUtils package.
# check which US-specific models are contained in the forecast archive.
model_names <- get_all_models(hub = "US")
There are 105 models within the Zoltar forecast archive from the U.S COVID-19 Forecasts Hub.
Load the incident cases forecasts models of 1 through 4 week horizons of a select number of US-specific models that are published by the CDC and are contained within the Zoltar forecast archive.
if (load_all_models == 1) {
# load forecasts of all models
system.time(all_model_inc_case_targets <- load_forecasts(
location = "US",
hub = "US",
types = c("point","quantile"),
targets = paste(1:4, "wk ahead inc case"),
as_of = "2021-09-06",
source = "zoltar"))
}
# load a select few point forecasts that were submitted in from zoltar forecast archive
inc_case_targets <- load_forecasts(
models = c("COVIDhub-ensemble", "CovidAnalytics-DELPHI", "CU-select",
"IHME-CurveFit", "LANL-GrowthRate","USC-SI_kJalpha",
"JHU_IDD-CovidSP", "UVA-Ensemble"),
location = "US",
hub = "US",
types = c("point","quantile"),
targets = paste(1:4, "wk ahead inc case"),
as_of = "2021-09-06",
source = "zoltar")
54 models stored in the Zoltar forecast archive have submitted a weekly incident case forecast.
# display top rows of forecast data
head(inc_case_targets)
inc_case_targets <- inc_case_targets %>%
mutate(forecast_day = lubridate::wday(forecast_date, label = TRUE, abbr = FALSE)) %>%
select(model, forecast_date, forecast_day, everything())
Filter the models for respective 1 through 4 week horizon point forecasts.
# filter for 1-4 week horizon point forecasts
for (wk in 1:4) {
wk_forecasts <- adj_inc_case_targets %>%
filter(horizon == wk,
type == "point") %>%
select(-quantile) %>%
dplyr::distinct(model, adj_forecast_date, .keep_all = TRUE)
assign(paste0("inc_case_targets_",wk,"_week"), wk_forecasts)
}
# inc_case_targets_1_week <- adj_inc_case_targets %>%
# filter(horizon == 1,
# type == "point") %>%
# select(-quantile)
Visualizing forecast data
Plot the CDC actual COVID-19 cases versus the forecast predictions for the select models.
# plot the weekly cases for United States according to CDC data as of 09/06/2021
p <- ggplot(data = national_cases_weekly, aes(x = week, y = new_cases)) +
geom_point() +
geom_line() +
labs(title = "CDC weekly incident COVID-19 cases")
p

# combine plots
p + geom_line(data = inc_case_targets_1_week, aes(x = target_end_date, y = value, color = model)) +
geom_point(data = inc_case_targets_1_week, aes(x = target_end_date, y = value, color = model)) +
labs(title = "CDC weekly incident COVID-19 cases versus various model point forecasts") +
theme(legend.position = "bottom")

Interactive plot of CDC observed cases versus 1 week horizon forecasts
fig1 <- plot_ly(type = 'scatter', mode = 'lines+markers')
fig1 <- fig1 %>%
add_trace(data = inc_case_targets_1_week, x = ~adj_forecast_date, y = ~value, color = ~factor(model)) %>%
add_trace(data = national_cases_weekly, x = ~week, y = ~new_cases, name = "CDC actual", color = I('black')) %>%
layout(title = "CDC weekly incident COVID-19 cases\n and various model point forecasts")
fig1
IHME-Curvefit only predicted incident COVID-cases for a few weeks before focusing on predicting hospitalizations and deaths, according to predictions stored in Zoltar forecast archive.
Identifying peaks of COVID-19 cases
# a 'peak' is defined as a local maxima with m points either side of it being smaller than it. hence, the bigger the parameter m, the more stringent is the peak finding procedure
# https://github.com/stas-g/findPeaks
find_peaks <- function (x, m = 3){
shape <- diff(sign(diff(x, na.pad = FALSE)))
pks <- sapply(which(shape < 0), FUN = function(i){
z <- i - m + 1
z <- ifelse(z > 0, z, 1)
w <- i + m + 1
w <- ifelse(w < length(x), w, length(x))
if(all(x[c(z : i, (i + 2) : w)] <= x[i + 1])) return(i + 1) else return(numeric(0))
})
pks <- unlist(pks)
pks
}
Peaks of observed weekly incident COVID-19 cases based on CDC data
# find peaks of observed cases
peak_cases <- national_cases_weekly[find_peaks(national_cases_weekly$new_cases, m = 3),]
# plot peaks for observed covid cases
q <- p + geom_point(data = peak_cases, aes(x = week, y = new_cases, color = "CDC actual")) +
labs(title = "Observed peaks of weekly incident COVID-19 cases",
x = "Weeks", y = "Incident cases", color = "Peaks") +
theme(legend.position = "bottom")
q

# get peak values of all values
model_peaks_list <- list()
for (i in unique(inc_case_targets_1_week$model)) {
models <- inc_case_targets_1_week %>%
filter(model == i)
peaks <- models[find_peaks(models$value, m = 3),]
model_peaks_list[[i]] <- peaks
}
# combine the df of peaks
model_peaks <- bind_rows(model_peaks_list)
Look at the peaks associated with forecast models
# combine plots of cdc actual peaks with forecast model peaks
q + geom_point(data = model_peaks, aes(x = adj_forecast_date, y = value, color = model))

# find magnitude and temporal differences in the peaks
cdc_peaks <- peak_cases %>%
mutate(model = rep("CDC", nrow(peak_cases))) %>%
select(week, model, new_cases) %>%
rename("peaks" = "new_cases")
model_observed_pks <- model_peaks %>%
select(adj_forecast_date, model, value) %>%
rename("week" = "adj_forecast_date",
"peaks" = "value") %>%
bind_rows(cdc_peaks)
Score each model weekly
To score models using covidHubUtils, data frame needs to be in same format as one gotten from load_truth function.
# load a truth data frame from covidHub
jhu_truth_df <- load_truth(
truth_source = c("JHU"),
target_variable = c("inc case"),
truth_end_date = "2021-09-04",
temporal_resolution = "weekly",
hub = "US",
locations = "US"
)
# top rows of truth dataframe
head(jhu_truth_df)
# to score models, data frame needs to be in same format as one gotten from load_truth
# TODO: convert national_cases_weekly to jhu_truth_df format
# TODO: use score_forecast function to compute weighted interval score
LS0tCnRpdGxlOiAiQ292aWQtMTkgZm9yZWNhc3QgbW9kZWwgZXhwbG9yYXRpb24iCnN1YnRpdGxlOiAiRXhwbG9yaW5nIGFuZCB2aXN1YWxpemluZyBDT1ZJRC0xOSBmb3JlY2FzdCBtb2RlbHMiCmRhdGU6ICJMYXN0IHVwZGF0ZWQ6IGByIGZvcm1hdChTeXMudGltZSgpLCclQiAlZCwgJVknKWAiCm91dHB1dDogCiAgaHRtbF9kb2N1bWVudDoKICAgIHRvYzogZmFsc2UKICAgIHRvY19mbG9hdDogZmFsc2UKICAgIGRmX3ByaW50OiBwYWdlZAogICAgY29kZV9kb3dubG9hZDogdHJ1ZQogICAgY29kZV9mb2xkaW5nOiBoaWRlCmVkaXRvcl9vcHRpb25zOiAKICBjaHVua19vdXRwdXRfdHlwZTogY29uc29sZQotLS0KCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQoZWNobyA9IFRSVUUsIGZpZy5hbGlnbiA9ICJjZW50ZXIiLCBmaWcud2lkdGg9MTIsIGZpZy5oZWlnaHQ9NikKYGBgCgpgYGB7ciwgZWNobz1GQUxTRSwgaW5jbHVkZT1GQUxTRX0KIyBjbGVhciBlbnZpcm9ubWVudApybShsaXN0ID0gbHMoKSkgICMgY2xlYXIgbWVtb3J5CmBgYAoKTG9hZCByZXF1aXJlZCBwYWNrYWdlcwpgYGB7ciwgbWVzc2FnZT1GQUxTRX0KIyBsb2FkIHBhY21hbiBwYWNrYWdlIHRvIGxvYWQgb3IgaW5zdGFsbCBvdGhlciByZXF1cmllZCBsaWJyYXJpZXMKaWYgKCFyZXF1aXJlKCdwYWNtYW4nKSkgaW5zdGFsbC5wYWNrYWdlcygncGFjbWFuJyk7IGxpYnJhcnkocGFjbWFuKQoKIyBsb2FkIChpbnN0YWxsIGlmIHJlcXVpcmVkKSBwYWNrYWdlcyBmcm9tIENSQU4KcF9sb2FkKCJoZXJlIiwidGlkeXZlcnNlIiwibHVicmlkYXRlIiwiamFuaXRvciIsInBsb3RseSIsImRvUGFyYWxsZWwiKQoKIyBsb2FkL2luc3RhbGwgcGFja2FnZXMgZnJvbSBHaXRIdWIKcF9sb2FkX2doKCJyZWljaGxhYi96b2x0ciIsICJyZWljaGxhYi9jb3ZpZEh1YlV0aWxzIikKYGBgCgpVdGlsaXppbmcgYHpvbHRyYCBwYWNrYWdlIHRvIGVuYWJsZSBhY2Nlc3MgdG8gWm9sdGFyIEFQSSB0byB0aGUgZm9yZWNhc3QgYXJjaGl2ZSBhbmQgYGNvdmlkSHViVXRpbHNgIHBhY2thZ2UgdG8gdGhhdCBwcm92aWRlIGZ1bmN0aW9ucyB0byByZWFkLCBwbG90IGFuZCBzY29yZSBmb3JlY2FzdCBkYXRhLgoKXCAgCgojIyMgSW1wb3J0aW5nIG9ic2VydmVkIGRhdGEKV2Ugd2lsbCB1c2UgdGhlIENvdmlkLTE5IGNhc2VzIGFzIHVwZGF0ZWQgYnkgQ0RDLiBXZSB3aWxsIGRvd25sb2FkIGRhdGEgZnJvbSB0aGUgZGF0YSBUYWJsZSBmb3IgZGFpbHkgQ2FzZSBUcmVuZHMgZm9yIFRoZSBVbml0ZWQgU3RhdGVzIFtmb3VuZCBoZXJlXSggaHR0cHM6Ly9jb3ZpZC5jZGMuZ292L2NvdmlkLWRhdGEtdHJhY2tlci8jdHJlbmRzX2RhaWx5Y2FzZXMpLiAKCkFsc28sIHdlIHdpbGwgcHVsbCBkYWlseSBDT1ZJRC0xOSBjYXNlcyBmcm9tIENEQyBhdCBhIHN0YXRlLWxldmVsIFtmb3VuZCBoZXJlXShodHRwczovL2RhdGEuY2RjLmdvdi9DYXNlLVN1cnZlaWxsYW5jZS9Vbml0ZWQtU3RhdGVzLUNPVklELTE5LUNhc2VzLWFuZC1EZWF0aHMtYnktU3RhdGUtby85bWZxLWNiMzYvZGF0YSkuCgpEYXRhIGV4dHJhY3RlZCBTZXB0ZW1iZXIsIDA2LCAyMDIxLiBIZXJlIGlzIHRoZSB0b3Agcm93cyBvZiB0aGUgaW5jaWRlbnQgQ09WSUQtMTkgY2FzZXMgZGF0YS4KYGBge3IsIG1lc3NhZ2U9RkFMU0V9CiMgY3N2IGV4cG9ydCBvZiBkYWlseSBjYXNlcyBmcm9tIENEQyBDT1ZJRCBkYXRhIHRyYWNrZXIgYXQgbmF0aW9uYWwgbGV2ZWwKbmF0aW9uYWxfY2FzZXNfZGFpbHkgPC0gcmVhZF9jc3YoaGVyZSgiZGF0YSIsICJkYXRhX3RhYmxlX2Zvcl9kYWlseV9jYXNlX3RyZW5kc19fdGhlX3VuaXRlZF9zdGF0ZXMuY3N2IiksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBjb2xfdHlwZXMgPSBjb2xzKERhdGUgPSBjb2xfZGF0ZShmb3JtYXQgPSAiJWIgJWQgJVkiKSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHNraXAgPSAyKQoKCiMgY3N2IGV4cG9ydCBvZiBkYWlseSBjYXNlcyBmcm9tIENEQyBDT1ZJRCBkYXRhIHRyYWNrZXIgYXQgc3RhdGUgbGV2ZWwKc3RhdGVzX2Nhc2VzX2RhaWx5IDwtIHJlYWRfY3N2KGhlcmUoImRhdGEiLCAiVW5pdGVkX1N0YXRlc19DT1ZJRC0xOV9DYXNlc19hbmRfRGVhdGhzX2J5X1N0YXRlX292ZXJfVGltZS5jc3YiKSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBjb2xfdHlwZXMgPSBjb2xzKHN1Ym1pc3Npb25fZGF0ZSA9IGNvbF9kYXRlKGZvcm1hdCA9ICIlbS8lZC8lWSIpKSkKCmhlYWQobmF0aW9uYWxfY2FzZXNfZGFpbHkpCmBgYAoKXCAgCgojIyMgUHJlcGFyaW5nIGRhdGEKCldlJ2xsIGFnZ3JlZ2F0ZSB0aGUgZGFpbHkgY2FzZSBjb3VudHMgaW50byB3ZWVrcywgc3RhcnRpbmcgb24gU3VuZGF5IGFzIGlzIHRoZSBzYW1lIGFzIGFuIGVwaWRlbWlvbG9naWNhbCB3ZWVrIChyZWZlcnJlZCB0byBhcyBNTVdSIHdlZWssIG1vcmUgaW5mbyBbZm91bmQgaGVyZV0oaHR0cHM6Ly9uZGMuc2VydmljZXMuY2RjLmdvdi93cC1jb250ZW50L3VwbG9hZHMvTU1XUl93ZWVrX292ZXJ2aWV3LnBkZikpCmBgYHtyfQojIGNsZWFuIHVwIGNvbHVtbiBuYW1lcwpuYXRpb25hbF9jYXNlc19kYWlseSA8LSBuYXRpb25hbF9jYXNlc19kYWlseSAlPiUKICBjbGVhbl9uYW1lcygpCgojIGFnZ3JlZ2F0ZSBkYWlseSBDT1ZJRCBjYXNlcyBpbnRvIHdlZWtseSBjYXNlIGNvdW50cywgCiMgc3RhcnQgd2VlayBvbiBTdW5kYXkgYWNjb3JkaW5nIHRvIGVwaWRlbWlvbG9naWNhbCB3ZWVrIChNTVdSIHdlZWspIAp3ZWVrIDwtIGFzLkRhdGUoY3V0KG5hdGlvbmFsX2Nhc2VzX2RhaWx5JGRhdGUsICJ3ZWVrIiwgc3RhcnQub24ubW9uZGF5ID0gRkFMU0UpKQpuYXRpb25hbF9jYXNlc193ZWVrbHkgPC0gYWdncmVnYXRlKG5ld19jYXNlcyB+IHdlZWssIG5hdGlvbmFsX2Nhc2VzX2RhaWx5LCBzdW0pCmhlYWQobmF0aW9uYWxfY2FzZXNfd2Vla2x5KQoKIyBmaWx0ZXIgZGF0ZXMgYWZ0ZXIgQXVndXN0Cm5hdGlvbmFsX2Nhc2VzX3dlZWtseSA8LSBuYXRpb25hbF9jYXNlc193ZWVrbHkgJT4lCiAgZmlsdGVyKHdlZWsgPCAiMjAyMS0wOS0wMSIpCmBgYAoKVGhlIGZvcmVjYXN0IG1vZGVscyBhcmUgcHVsbGVkIGZyb20gdGhlIFtab2x0YXIgZm9yZWNhc3QgbW9kZWwgYXJjaGl2ZV0oaHR0cHM6Ly93d3cuem9sdGFyZGF0YS5jb20vKSB1c2luZyB0aGUgW2NvdmlkSHViVXRpbHMgcGFja2FnZV0oaHR0cDovL3JlaWNobGFiLmlvL2NvdmlkSHViVXRpbHMvaW5kZXguaHRtbCkuIAoKYGBge3IsIG1lc3NhZ2U9RkFMU0V9CiMgY2hlY2sgd2hpY2ggVVMtc3BlY2lmaWMgbW9kZWxzIGFyZSBjb250YWluZWQgaW4gdGhlIGZvcmVjYXN0IGFyY2hpdmUuCm1vZGVsX25hbWVzIDwtIGdldF9hbGxfbW9kZWxzKGh1YiA9ICJVUyIpIApgYGAKClRoZXJlIGFyZSBgciBsZW5ndGgobW9kZWxfbmFtZXMpYCBtb2RlbHMgd2l0aGluIHRoZSBab2x0YXIgZm9yZWNhc3QgYXJjaGl2ZSBmcm9tIHRoZSBVLlMgQ09WSUQtMTkgRm9yZWNhc3RzIEh1Yi4KCkxvYWQgdGhlIGluY2lkZW50IGNhc2VzIGZvcmVjYXN0cyBtb2RlbHMgb2YgMSB0aHJvdWdoIDQgd2VlayBob3Jpem9ucyBvZiBhIHNlbGVjdCBudW1iZXIgb2YgVVMtc3BlY2lmaWMgbW9kZWxzIHRoYXQgYXJlIHB1Ymxpc2hlZCBieSB0aGUgQ0RDIGFuZCBhcmUgY29udGFpbmVkIHdpdGhpbiB0aGUgWm9sdGFyIGZvcmVjYXN0IGFyY2hpdmUuCgpgYGB7ciwgaW5jbHVkZT1GQUxTRX0KIyBpbmNsdWRlIGFsbCBmb3JlY2FzdCBtb2RlbHMgb3IgYSBmZXcgc3BlY2lmaWVkIG1vZGVscwpsb2FkX2FsbF9tb2RlbHMgPC0gMApgYGAKCmBgYHtyLCBtZXNzYWdlPUZBTFNFLCBjYWNoZT1UUlVFfQppZiAobG9hZF9hbGxfbW9kZWxzID09IDEpIHsKIyBsb2FkIGZvcmVjYXN0cyBvZiBhbGwgbW9kZWxzICAKc3lzdGVtLnRpbWUoYWxsX21vZGVsX2luY19jYXNlX3RhcmdldHMgPC0gbG9hZF9mb3JlY2FzdHMoCiAgbG9jYXRpb24gPSAiVVMiLAogIGh1YiA9ICJVUyIsCiAgdHlwZXMgPSBjKCJwb2ludCIsInF1YW50aWxlIiksCiAgdGFyZ2V0cyA9ICBwYXN0ZSgxOjQsICJ3ayBhaGVhZCBpbmMgY2FzZSIpLAogIGFzX29mID0gIjIwMjEtMDktMDYiLAogIHNvdXJjZSA9ICJ6b2x0YXIiKSkKICB9CmBgYAoKYGBge3IsIG1lc3NhZ2U9RkFMU0UsIGNhY2hlPVRSVUV9CiMgbG9hZCBhIHNlbGVjdCBmZXcgcG9pbnQgZm9yZWNhc3RzIHRoYXQgd2VyZSBzdWJtaXR0ZWQgaW4gZnJvbSB6b2x0YXIgZm9yZWNhc3QgYXJjaGl2ZQogIGluY19jYXNlX3RhcmdldHMgPC0gbG9hZF9mb3JlY2FzdHMoCiAgbW9kZWxzID0gYygiQ09WSURodWItZW5zZW1ibGUiLCAiQ292aWRBbmFseXRpY3MtREVMUEhJIiwgIkNVLXNlbGVjdCIsCiAgICAgICAgICAgICAiSUhNRS1DdXJ2ZUZpdCIsICJMQU5MLUdyb3d0aFJhdGUiLCJVU0MtU0lfa0phbHBoYSIsIAogICAgICAgICAgICAgIkpIVV9JREQtQ292aWRTUCIsICJVVkEtRW5zZW1ibGUiKSwKICBsb2NhdGlvbiA9ICJVUyIsCiAgaHViID0gIlVTIiwKICB0eXBlcyA9IGMoInBvaW50IiwicXVhbnRpbGUiKSwKICB0YXJnZXRzID0gIHBhc3RlKDE6NCwgIndrIGFoZWFkIGluYyBjYXNlIiksCiAgYXNfb2YgPSAiMjAyMS0wOS0wNiIsCiAgc291cmNlID0gInpvbHRhciIpCmBgYAoKYGBge3IsIGluY2x1ZGU9RkFMU0V9CiMgdGhlIG51bWJlciBvZiBtb2RlbHMgd2l0aCB3ZWVrIGluY2lkZW50IGNhc2UgZm9yZWNhc3RzCmlmIChsb2FkX2FsbF9tb2RlbHMgPT0gMSl7CiAgbGVuZ3RoKHVuaXF1ZShhbGxfbW9kZWxfaW5jX2Nhc2VfdGFyZ2V0cyRtb2RlbCkpCiAgfSBlbHNlIHsKICAgIGxlbmd0aCh1bmlxdWUoaW5jX2Nhc2VfdGFyZ2V0cyRtb2RlbCkpCiAgfQpgYGAKCjU0IG1vZGVscyBzdG9yZWQgaW4gdGhlIFpvbHRhciBmb3JlY2FzdCBhcmNoaXZlIGhhdmUgc3VibWl0dGVkIGEgd2Vla2x5IGluY2lkZW50IGNhc2UgZm9yZWNhc3QuCgpgYGB7cn0KIyBkaXNwbGF5IHRvcCByb3dzIG9mIGZvcmVjYXN0IGRhdGEKaGVhZChpbmNfY2FzZV90YXJnZXRzKQoKaW5jX2Nhc2VfdGFyZ2V0cyA8LSBpbmNfY2FzZV90YXJnZXRzICU+JQogIG11dGF0ZShmb3JlY2FzdF9kYXkgPSBsdWJyaWRhdGU6OndkYXkoZm9yZWNhc3RfZGF0ZSwgbGFiZWwgPSBUUlVFLCBhYmJyID0gRkFMU0UpKSAlPiUKICBzZWxlY3QobW9kZWwsIGZvcmVjYXN0X2RhdGUsIGZvcmVjYXN0X2RheSwgZXZlcnl0aGluZygpKQpgYGAKCmBgYHtyLCBpbmNsdWRlPUZBTFNFfQojIGFkanVzdCBmb3JlY2FzdCBkYXRlcyB0byBsaWUgb24gdGhlIHNhbWUgd2VlayBvZiB0aGVpciByZXNwZWN0aXZlIGZvcmVjYXN0IHRhcmdldCBlbmQgZGF0ZQojIGh0dHBzOi8vZ2l0aHViLmNvbS9yZWljaGxhYi9jb3ZpZDE5LWZvcmVjYXN0LWh1Yi9ibG9iL21hc3Rlci9kYXRhLXByb2Nlc3NlZC9SRUFETUUubWQjZm9yZWNhc3QtZmlsZS1mb3JtYXQKYWRqX2luY19jYXNlX3RhcmdldHMgPC0gaW5jX2Nhc2VfdGFyZ2V0cyAlPiUKICBtdXRhdGUoYWRqX2ZvcmVjYXN0X2RhdGUgPSBhcy5EYXRlKGlmZWxzZShmb3JlY2FzdF9kYXkgPT0gIlN1bmRheSIsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmb3JlY2FzdF9kYXRlLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBjYXNlX3doZW4oZm9yZWNhc3RfZGF5ID09ICJNb25kYXkiIH4gZm9yZWNhc3RfZGF0ZSAtIDEsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmb3JlY2FzdF9kYXkgPT0gIlR1ZXNkYXkiIH4gZm9yZWNhc3RfZGF0ZSArIDUsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmb3JlY2FzdF9kYXkgPT0gIldlZG5lc2RheSIgfiBmb3JlY2FzdF9kYXRlICsgNCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZvcmVjYXN0X2RheSA9PSAiVGh1cnNkYXkiIH4gZm9yZWNhc3RfZGF0ZSArIDMsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmb3JlY2FzdF9kYXkgPT0gIkZyaWRheSIgfiBmb3JlY2FzdF9kYXRlICsgMiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZvcmVjYXN0X2RheSA9PSAiU2F0dXJkYXkiIH4gZm9yZWNhc3RfZGF0ZSArIDEpKSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG9yaWdpbiA9ICIxOTcwLTAxLTAxIiksCiAgICAgICAgIGFkal9mb3JlY2FzdF9kYXkgPSBsdWJyaWRhdGU6OndkYXkoYWRqX2ZvcmVjYXN0X2RhdGUsIGxhYmVsID0gVFJVRSwgYWJiciA9IEZBTFNFKSkgJT4lCiAgc2VsZWN0KGZvcmVjYXN0X2RhdGUsIGZvcmVjYXN0X2RheSwgYWRqX2ZvcmVjYXN0X2RhdGUsIGFkal9mb3JlY2FzdF9kYXksIHRhcmdldF9lbmRfZGF0ZSwgZXZlcnl0aGluZygpKSAlPiUKICBhcnJhbmdlKGFkal9mb3JlY2FzdF9kYXRlKSAlPiUKICAjIGRwbHlyOjpkaXN0aW5jdChtb2RlbCwgYWRqX2ZvcmVjYXN0X2RhdGUsIC5rZWVwX2FsbCA9IFRSVUUpICU+JQogIGZpbHRlcihhZGpfZm9yZWNhc3RfZGF0ZSA8ICIyMDIxLTA5LTAxIikKYGBgCgpGaWx0ZXIgdGhlIG1vZGVscyBmb3IgcmVzcGVjdGl2ZSAxIHRocm91Z2ggNCB3ZWVrIGhvcml6b24gcG9pbnQgZm9yZWNhc3RzLgpgYGB7cn0KIyBmaWx0ZXIgZm9yIDEtNCB3ZWVrIGhvcml6b24gcG9pbnQgZm9yZWNhc3RzCmZvciAod2sgaW4gMTo0KSB7CiAgd2tfZm9yZWNhc3RzIDwtIGFkal9pbmNfY2FzZV90YXJnZXRzICU+JSAKICBmaWx0ZXIoaG9yaXpvbiA9PSB3aywKICAgICAgICAgdHlwZSA9PSAicG9pbnQiKSAlPiUKICBzZWxlY3QoLXF1YW50aWxlKSAlPiUgCiAgZHBseXI6OmRpc3RpbmN0KG1vZGVsLCBhZGpfZm9yZWNhc3RfZGF0ZSwgLmtlZXBfYWxsID0gVFJVRSkKICBhc3NpZ24ocGFzdGUwKCJpbmNfY2FzZV90YXJnZXRzXyIsd2ssIl93ZWVrIiksIHdrX2ZvcmVjYXN0cykKfQoKIyBpbmNfY2FzZV90YXJnZXRzXzFfd2VlayA8LSBhZGpfaW5jX2Nhc2VfdGFyZ2V0cyAlPiUgCiMgICBmaWx0ZXIoaG9yaXpvbiA9PSAxLAojICAgICAgICAgIHR5cGUgPT0gInBvaW50IikgJT4lCiMgICBzZWxlY3QoLXF1YW50aWxlKQpgYGAKCiMjIyBWaXN1YWxpemluZyBmb3JlY2FzdCBkYXRhCgpQbG90IHRoZSBDREMgYWN0dWFsIENPVklELTE5IGNhc2VzIHZlcnN1cyB0aGUgZm9yZWNhc3QgcHJlZGljdGlvbnMgZm9yIHRoZSBzZWxlY3QgbW9kZWxzLgpgYGB7cn0KIyBwbG90IHRoZSB3ZWVrbHkgY2FzZXMgZm9yIFVuaXRlZCBTdGF0ZXMgYWNjb3JkaW5nIHRvIENEQyBkYXRhIGFzIG9mIDA5LzA2LzIwMjEKcCA8LSBnZ3Bsb3QoZGF0YSA9IG5hdGlvbmFsX2Nhc2VzX3dlZWtseSwgYWVzKHggPSB3ZWVrLCB5ID0gbmV3X2Nhc2VzKSkgKwogIGdlb21fcG9pbnQoKSArCiAgZ2VvbV9saW5lKCkgKwogIGxhYnModGl0bGUgPSAiQ0RDIHdlZWtseSBpbmNpZGVudCBDT1ZJRC0xOSBjYXNlcyIpCnAKYGBgCgpgYGB7ciwgaW5jbHVkZT1GQUxTRX0KIyBwbG90IHZhcmlvdXMgbW9kZWwgMSB3ZWVrIGhvcml6b24gZm9yZWNhc3RzCmdncGxvdChkYXRhID0gaW5jX2Nhc2VfdGFyZ2V0c18xX3dlZWssIGFlcyh4ID0gYWRqX2ZvcmVjYXN0X2RhdGUsIHkgPSB2YWx1ZSwgY29sb3IgPSBtb2RlbCkpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21fbGluZSgpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAiYm90dG9tIikKYGBgCgpgYGB7cn0KIyBjb21iaW5lIHBsb3RzCnAgKyBnZW9tX2xpbmUoZGF0YSA9IGluY19jYXNlX3RhcmdldHNfMV93ZWVrLCBhZXMoeCA9IHRhcmdldF9lbmRfZGF0ZSwgeSA9IHZhbHVlLCBjb2xvciA9IG1vZGVsKSkgKwogIGdlb21fcG9pbnQoZGF0YSA9IGluY19jYXNlX3RhcmdldHNfMV93ZWVrLCBhZXMoeCA9IHRhcmdldF9lbmRfZGF0ZSwgeSA9IHZhbHVlLCBjb2xvciA9IG1vZGVsKSkgKwogIGxhYnModGl0bGUgPSAiQ0RDIHdlZWtseSBpbmNpZGVudCBDT1ZJRC0xOSBjYXNlcyB2ZXJzdXMgdmFyaW91cyBtb2RlbCBwb2ludCBmb3JlY2FzdHMiKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gImJvdHRvbSIpCmBgYAoKSW50ZXJhY3RpdmUgcGxvdCBvZiBDREMgb2JzZXJ2ZWQgY2FzZXMgdmVyc3VzIDEgd2VlayBob3Jpem9uIGZvcmVjYXN0cwpgYGB7cn0KZmlnMSA8LSBwbG90X2x5KHR5cGUgPSAnc2NhdHRlcicsICBtb2RlID0gJ2xpbmVzK21hcmtlcnMnKQpmaWcxIDwtIGZpZzEgJT4lIAogIGFkZF90cmFjZShkYXRhID0gaW5jX2Nhc2VfdGFyZ2V0c18xX3dlZWssIHggPSB+YWRqX2ZvcmVjYXN0X2RhdGUsIHkgPSB+dmFsdWUsIGNvbG9yID0gfmZhY3Rvcihtb2RlbCkpICU+JQogIGFkZF90cmFjZShkYXRhID0gbmF0aW9uYWxfY2FzZXNfd2Vla2x5LCB4ID0gfndlZWssIHkgPSB+bmV3X2Nhc2VzLCBuYW1lID0gIkNEQyBhY3R1YWwiLCBjb2xvciA9IEkoJ2JsYWNrJykpICU+JQogIGxheW91dCh0aXRsZSA9ICJDREMgd2Vla2x5IGluY2lkZW50IENPVklELTE5IGNhc2VzXG4gYW5kIHZhcmlvdXMgbW9kZWwgcG9pbnQgZm9yZWNhc3RzIikKCmZpZzEKYGBgCgpJSE1FLUN1cnZlZml0IG9ubHkgcHJlZGljdGVkIGluY2lkZW50IENPVklELWNhc2VzIGZvciBhIGZldyB3ZWVrcyBiZWZvcmUgZm9jdXNpbmcgb24gcHJlZGljdGluZyBob3NwaXRhbGl6YXRpb25zIGFuZCBkZWF0aHMsIGFjY29yZGluZyB0byBwcmVkaWN0aW9ucyBzdG9yZWQgaW4gWm9sdGFyIGZvcmVjYXN0IGFyY2hpdmUuCgpJZGVudGlmeWluZyBwZWFrcyBvZiBDT1ZJRC0xOSBjYXNlcwpgYGB7cn0KIyBhICdwZWFrJyBpcyBkZWZpbmVkIGFzIGEgbG9jYWwgbWF4aW1hIHdpdGggbSBwb2ludHMgZWl0aGVyIHNpZGUgb2YgaXQgYmVpbmcgc21hbGxlciB0aGFuIGl0LiBoZW5jZSwgdGhlIGJpZ2dlciB0aGUgcGFyYW1ldGVyIG0sIHRoZSBtb3JlIHN0cmluZ2VudCBpcyB0aGUgcGVhayBmaW5kaW5nIHByb2NlZHVyZQojIGh0dHBzOi8vZ2l0aHViLmNvbS9zdGFzLWcvZmluZFBlYWtzCmZpbmRfcGVha3MgPC0gZnVuY3Rpb24gKHgsIG0gPSAzKXsKICAgIHNoYXBlIDwtIGRpZmYoc2lnbihkaWZmKHgsIG5hLnBhZCA9IEZBTFNFKSkpCiAgICBwa3MgPC0gc2FwcGx5KHdoaWNoKHNoYXBlIDwgMCksIEZVTiA9IGZ1bmN0aW9uKGkpewogICAgICAgeiA8LSBpIC0gbSArIDEKICAgICAgIHogPC0gaWZlbHNlKHogPiAwLCB6LCAxKQogICAgICAgdyA8LSBpICsgbSArIDEKICAgICAgIHcgPC0gaWZlbHNlKHcgPCBsZW5ndGgoeCksIHcsIGxlbmd0aCh4KSkKICAgICAgIGlmKGFsbCh4W2MoeiA6IGksIChpICsgMikgOiB3KV0gPD0geFtpICsgMV0pKSByZXR1cm4oaSArIDEpIGVsc2UgcmV0dXJuKG51bWVyaWMoMCkpCiAgICB9KQogICAgIHBrcyA8LSB1bmxpc3QocGtzKQogICAgIHBrcwp9CmBgYAoKClBlYWtzIG9mIG9ic2VydmVkIHdlZWtseSBpbmNpZGVudCBDT1ZJRC0xOSBjYXNlcyBiYXNlZCBvbiBDREMgZGF0YQpgYGB7cn0KIyBmaW5kIHBlYWtzIG9mIG9ic2VydmVkIGNhc2VzCnBlYWtfY2FzZXMgPC0gbmF0aW9uYWxfY2FzZXNfd2Vla2x5W2ZpbmRfcGVha3MobmF0aW9uYWxfY2FzZXNfd2Vla2x5JG5ld19jYXNlcywgbSA9IDMpLF0KCiMgcGxvdCBwZWFrcyBmb3Igb2JzZXJ2ZWQgY292aWQgY2FzZXMKcSA8LSBwICsgZ2VvbV9wb2ludChkYXRhID0gcGVha19jYXNlcywgYWVzKHggPSB3ZWVrLCB5ID0gbmV3X2Nhc2VzLCBjb2xvciA9ICJDREMgYWN0dWFsIikpICsKICBsYWJzKHRpdGxlID0gIk9ic2VydmVkIHBlYWtzIG9mIHdlZWtseSBpbmNpZGVudCBDT1ZJRC0xOSBjYXNlcyIsIAogICAgICAgeCA9ICJXZWVrcyIsIHkgPSAiSW5jaWRlbnQgY2FzZXMiLCBjb2xvciA9ICJQZWFrcyIpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAiYm90dG9tIikKCnEKYGBgCgpgYGB7cn0KIyBnZXQgcGVhayB2YWx1ZXMgb2YgYWxsIHZhbHVlcwptb2RlbF9wZWFrc19saXN0IDwtIGxpc3QoKQpmb3IgKGkgaW4gdW5pcXVlKGluY19jYXNlX3RhcmdldHNfMV93ZWVrJG1vZGVsKSkgewogIG1vZGVscyA8LSBpbmNfY2FzZV90YXJnZXRzXzFfd2VlayAlPiUKICAgIGZpbHRlcihtb2RlbCA9PSBpKQogIHBlYWtzIDwtIG1vZGVsc1tmaW5kX3BlYWtzKG1vZGVscyR2YWx1ZSwgbSA9IDMpLF0KICBtb2RlbF9wZWFrc19saXN0W1tpXV0gPC0gcGVha3MKfQoKIyBjb21iaW5lIHRoZSBkZiBvZiBwZWFrcwptb2RlbF9wZWFrcyA8LSBiaW5kX3Jvd3MobW9kZWxfcGVha3NfbGlzdCkKYGBgCgpMb29rIGF0IHRoZSBwZWFrcyBhc3NvY2lhdGVkIHdpdGggZm9yZWNhc3QgbW9kZWxzCmBgYHtyfQojIGNvbWJpbmUgcGxvdHMgb2YgY2RjIGFjdHVhbCBwZWFrcyB3aXRoIGZvcmVjYXN0IG1vZGVsIHBlYWtzCnEgKyBnZW9tX3BvaW50KGRhdGEgPSBtb2RlbF9wZWFrcywgYWVzKHggPSBhZGpfZm9yZWNhc3RfZGF0ZSwgeSA9IHZhbHVlLCBjb2xvciA9IG1vZGVsKSkKYGBgCgpgYGB7cn0KIyBmaW5kIG1hZ25pdHVkZSBhbmQgdGVtcG9yYWwgZGlmZmVyZW5jZXMgaW4gdGhlIHBlYWtzCmNkY19wZWFrcyA8LSBwZWFrX2Nhc2VzICU+JQogIG11dGF0ZShtb2RlbCA9IHJlcCgiQ0RDIiwgbnJvdyhwZWFrX2Nhc2VzKSkpICU+JQogIHNlbGVjdCh3ZWVrLCBtb2RlbCwgbmV3X2Nhc2VzKSAlPiUKICByZW5hbWUoInBlYWtzIiA9ICJuZXdfY2FzZXMiKQoKbW9kZWxfb2JzZXJ2ZWRfcGtzIDwtIG1vZGVsX3BlYWtzICU+JQogIHNlbGVjdChhZGpfZm9yZWNhc3RfZGF0ZSwgbW9kZWwsIHZhbHVlKSAlPiUKICByZW5hbWUoIndlZWsiID0gImFkal9mb3JlY2FzdF9kYXRlIiwKICAgICAgICAgInBlYWtzIiA9ICJ2YWx1ZSIpICU+JQogIGJpbmRfcm93cyhjZGNfcGVha3MpCgpgYGAKCmBgYHtyLCBpbmNsdWRlPUZBTFNFfQojIGZpbmQgZWFybGllc3QgcGVhayBmb3IgZWFjaCBvZiB0aGUgbW9kZWxzCm1vZGVsX29ic2VydmVkX3BrcyAlPiUKICBncm91cF9ieShtb2RlbCkgJT4lCiAgc3VtbWFyaXNlKGZpcnN0X3BrX2RhdGUgPSBtaW4od2VlaykpCmBgYAoKIyMjIFNjb3JlIGVhY2ggbW9kZWwgd2Vla2x5ClRvIHNjb3JlIG1vZGVscyB1c2luZyBjb3ZpZEh1YlV0aWxzLCBkYXRhIGZyYW1lIG5lZWRzIHRvIGJlIGluIHNhbWUgZm9ybWF0IGFzIG9uZSBnb3R0ZW4gZnJvbSBgbG9hZF90cnV0aGAgZnVuY3Rpb24uCmBgYHtyLCBtZXNzYWdlPUZBTFNFfQojIGxvYWQgYSB0cnV0aCBkYXRhIGZyYW1lIGZyb20gY292aWRIdWIKamh1X3RydXRoX2RmIDwtIGxvYWRfdHJ1dGgoCiAgdHJ1dGhfc291cmNlID0gYygiSkhVIiksCiAgdGFyZ2V0X3ZhcmlhYmxlID0gYygiaW5jIGNhc2UiKSwKICB0cnV0aF9lbmRfZGF0ZSA9ICIyMDIxLTA5LTA0IiwKICB0ZW1wb3JhbF9yZXNvbHV0aW9uID0gIndlZWtseSIsCiAgaHViID0gIlVTIiwKICBsb2NhdGlvbnMgPSAgIlVTIgopCgojIHRvcCByb3dzIG9mIHRydXRoIGRhdGFmcmFtZQpoZWFkKGpodV90cnV0aF9kZikKYGBgCgpgYGB7cn0KIyB0byBzY29yZSBtb2RlbHMsIGRhdGEgZnJhbWUgbmVlZHMgdG8gYmUgaW4gc2FtZSBmb3JtYXQgYXMgb25lIGdvdHRlbiBmcm9tIGxvYWRfdHJ1dGgKIyBUT0RPOiBjb252ZXJ0IG5hdGlvbmFsX2Nhc2VzX3dlZWtseSB0byBqaHVfdHJ1dGhfZGYgZm9ybWF0CiMgVE9ETzogdXNlIHNjb3JlX2ZvcmVjYXN0IGZ1bmN0aW9uIHRvIGNvbXB1dGUgd2VpZ2h0ZWQgaW50ZXJ2YWwgc2NvcmUKYGBgCgoK